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Abstract 

In spiking neural networks, the information is conveyed by the spike times, that depend on the intrinsic 
dynamics of each neuron, the input they receive and on the connections between neurons. In this article 
we study the Markovian nature of the sequence of spike times in stochastic neural networks, and in 
particular the ability to deduce from a spike train the next spike time, and therefore produce a description 
of the network activity only based on the spike times regardless of the membrane potential process. 

To study this question in a rigorous manner, we introduce and study an event-based description of 
networks of noisy integrate-and-fire neurons, i.e. that is based on the computation of the spike times. 
We show that the firing times of the neurons in the networks constitute a Markov chain, whose transition 
probability is related to the probability distribution of the interspike interval of the neurons in the network. 
In the cases where the Markovian model can be developed, the transition probability is explicitly derived 
in such classical cases of neural networks as the linear integrate-and-fire neuron models with excitatory 
and inhibitory interactions, for different types of synapses, possibly featuring noisy synaptic integration, 
transmission delays and absolute and relative refractory period. This covers most of the cases that have 
been investigated in the event-based description of spiking deterministic neural networks. 



Introduction 



Growing experimental evi dence tends to establish tha t 



ral computations (see e.g. ( IDelorme and Thorpe , 



2001 



spike timings are essential to account for ne u 



Fabre-Thorpe et al 



1998 



Thorpe et al 



2001)). 



This fact has motivated the use of spiking neuron models, rather than the traditional rate-based models, 
and a very general interest for the characterization of the structure of spike times as well as for their 
precise and efficient computation. 

One of the great challenges in computational neuroscience consists in inferring from a sequence of 
spike times the next spike that will be fired in the network. In other words, knowing the input a neuron 



1 



receives and the network's connectivity map, can one compute the dynamics of the spike times regardless 
of the actual value of the underlying membrane potential value? Answering this qu estion amounts defi n- 



2007)). 



ing what is called an event-based description of the network (see e.g. the review of dBrette et al. 
The purpose of these models is to provide an exact simulation of the network, and this is precisely based 
on the ability to infer the next spike time based on the knowledge of the previous spike times and on the 
parameters of the m odels, without refer ring to the underlying membrane potential. 



As discussed in dBrette et al 



20071) . in the deterministic case, this property of event-based models is 



true only for very few cases of simple types of neuron models, where explicit solutions to the membrane 
potential voltage equation are available. Deter ministic even t -based algorithms were only developed for 



simpl e pulse-coupled integrate-and-fire models dWattsL 



1994 



Claverol et al 



2002; 



2002; 



2003 ) , and more complex ones in the case of the Spike Response Model dMakin 



Gerstner and Kistler 



to. 



Delorme and Thorpe, 



2003 



Marian et al 



2002b ). Rudo l ph and Destexhe developed a class of models that allow event- 



based simulation (Rudolph and Destexhe, 



Brette described in dBrette , 



2006, 



20061) including a model of synaptic conductances. Romain 



20071) algorithms that apply to some classes of input current (exponen- 



tial currents), and Tonnelier and collaborators extended these works to the c ase of the quadratic int egrate 



and-fire models in a case where the equation is integrable in closed-form dTonnelier et al 



2007). Even 



in this case the authors must make quite restrictive assumptions with regards to the connectivity and the 
input current while being more general in terms of the intrinsic dynamics. Their approach relies on the 
fact that the solution of the particular equation they choose has a closed-form solution. The price to pay is 
that no generalization is possible. All these models share the property that the solution to the membrane- 
potential equation can be found in closed-form, and therefore the spike times can be expressed as the 
times of the first intersection of these solutions to the threshold, i.e. by inverting explicitly the solution. 
Moreover, these models only apply to deterministic inputs of a very specific form (either constant, or 
some types of exponential currents for the linear IF neuron). 

They also have the common property that they do not take into account the natural randomness present 
in the input to the neurons. In the cortex, noise constitutes a prom inent aspect of the neur al activity, and 
is even sometimes assumed to have a functional role (see e.g. dRolls and Decol 120101) ). A common 
model for this random arrival of spikes consists in considering that the neuron receives a Brownian noisy 
current (diffusion approximation), as we further explain in section [T] Taking into account the effect 
of noise in the spike timings is a very difficult task, even for single neurons, and has been ex tensively 
studie d in the past decades and sti ll con stitutes an active area of rese arch (see e.g. the books dH olden, 



19761) . (IRicciardi and Smith , 



19771) and dGerstner and Kistlei , 



2002bl) or dTouboul and Faugeras , 



2007) 



and the references therein for a review). 

In the stochastic case, event-based descriptions of spike timings become much more complicated. 
Indeed, even when we have a closed-form expression of the membrane potential (which only occurs in 
very few specific cases) there are at least two reasons why the problem is extremely difficult to solve. 
First, it is not always possible to invert the probability density function of the membrane potential to find 
the first crossing time to the threshold. Second is the fact that the spike times of all neurons can have a 
very complex correlation structure with the underlying membrane potential voltage: even if we were able 
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to compute the probability distribution of the next spike numerically, the dependence on the membrane 
potential variable makes the next spike time extremely difficult to tame. 

Extending these results to the statistics of spike times in a network constitutes therefore the next step 
in both the theory of event-based descriptions of networks and of spike statistics. Indeed, the approach 
proposed in the deterministic case of solving in closed-form the membrane-potential equatio n and invert 



ing it will fail. This problem was already identified in the review by Brette and collaborators (IBrette et al 



2007), where the authors mention that "the problem of simulating directly a stochastic process in asyn- 



chronous algorithms is much harder because even for the simplest stochastic neuron models, there is no 
closed analytical formula for the distribution of the time to the next spike (see e.g. ( lTuckwelllll988l) ')". 



Another extremely complex problem lies in the dependence structure of the spike times among the differ- 
ent neurons of the network, and the statistical correlation of the next spike time to the trajectory (sample 
path) of the membrane potential. These are precisely the problems we address in the present manuscript. 

Existing work on event-based descriptions of neuronal networks is based on a perfect match between 
spike times of the voltage process and the event-based model. If this aim is highly relevant in the deter- 
ministic case, in the stochastic case the relevant information is the probability distribution of the spike 
times (i.e. the probability of a given neuron to fire at a given time) rather than some specific realizations 
of spike trains. Our approach therefore aims at perfectly reproducing the probability distribution of the 
spike times of the network. 

In all the paper, we consider networks composed of TV interacting integrate-and-fire spiking neurons 
subjected to noisy inputs from the external environment. Mathematically, the membrane potential of 
these neurons is defined by a stochastic process: each neuron receives at its synapses noisy inputs corre- 
sponding to the random activity of ion channels and to the external activity of the network. During the 
time intervals where no spike is emitted in the network, the membrane potential of each neuron evolves 
independently according to its intrinsic dynamics. When the membrane potential (t) of the neuron 
indexed by i reaches its deterministic threshold function 6 if) at time to, the neuron elicits an action 
potential. Subsequently, its membrane potential is reset to a given value V r , and the states of all its 
postsynaptic neurons j are modified. We denote by V(i) the postsynaptic neighborhood of the neuron i, 
i.e. the set of neurons that receive the spikes fired by neuron i. When a spike is transmitted to a neigh- 
boring cell, it results in the emergence of a positive or negative postsynaptic pulse whose sign depends 
on the connection type and whose shape depends on the neurotransmitters and the receptors. 

In this article, we address the problem of event-based models for stochastic networks, and more 
fundamentally of the Markovian nature of the spike times, in some simple types of networks of integrate- 
and-fire neurons with different intrinsic dynamics and different kinds of synaptic integration. To this 
purpose, we introduce a special mathematical framework to describe the dynamics of the spikes in sec- 
tionQ] We then start by addressing the case of purely inhibitory networks, and we show in section [2] that 
the spike times can be described through the introduction of a Markov chain whose transition probability 
is related to a first passage-time problem, for some simple types of neuron models. We further extend 
this approach in appendix|A]to more complex bio-inspired neuronal descriptions, and show the Markov 
property does not necessarily hold true in all cases. The more slightly complex case of networks includ- 
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ing excitatory interactions is then treated in section [3] and the same conclusion holds true: a Markov 
chain can be defined in relation with the spike times in the network, whose transitions are related to first 
passage time problems. The present framework and derivations open the way to different theoretical, 
numerical and computational developments, some of which are discussed in the conclusion. 



1 Theoretical framework 

We start by introducing the general mathematical concepts and the framework that we will be using in this 
manuscript in order to study spike times in stochastic integrate-and-fire networks. The first and central 
concept we want to recall is the Markov property for a random sequence. Basically, assume that the 
random sequence (-Xfc)fceN is defined by a certain recurrence property given by the transition probability 
density p(X n+ i\X n , . . . , Xq). This expresses probability density of the process X n+ i conditionally to 
the values of Xq, ■ ■ ■ , X n . The random sequence is called a Markov chain if this conditional probability 
depends only on X n , being independent of (Xk)k=o,...n-i, that is: 

p(X n+ i = x n+1 \X n = x„, Xq = xq) = p(X n+ i = x n+ i\X n = x n ). 

Loosely speaking, the knowledge of the present state makes the future of the process independent of its 
past states. 

The building block of our framework consists in considering the dynamics of the spike times instead 
of the dynamics of all variables involved in the description of the neuron's state. This purpose motivates 
the introduction of what we call the countdown process, a stochastic process governing the dynamics of 
the spikes in the network and defined as follows: 

Definition 1.1. [Countdown process] For each neuron i, we define X^ (t) > to be the remaining time 
until the next emission of a spike by neuron i if it does not receive any spike meanwhile. We call this 
stochastic process the countdown process of the neuron i. Following the tradition established in the 
stochastic calculus literature, we also sometimes note x[ for X^ l \t). The N dimensional process 
X t = (X^\ . . . , X^ N) ) is called the countdown process of the network, or simply countdown process. 

For each neuron i G {1, . . ., N}, the dynamics of (XW(i)) is very simple, linearly decreasing with 
slope —1 during the intervals of time where no spike is received or produced: 

dX®(t) 

txjtL—i (1) 

At time t, the next spike will be fired by neuron i = Arg Min jel N X^> (i) at time t + X^ (t). 

At spike time, the membrane potential of the neuron that just fired an action potential is reset, and 
therefore the new countdown value after the spike emission is updated to a new value corresponding to 
the next spike time of this neuron if nothing occurs meanwhile. This random variable, noted yW, is 
referred to as the reset random variable. 

Each neuron connected to the one that fired the action potential has his state modified, and its next 
firing time is reset accordingly, conditionally to the value of the next firing time predicted before the 
reception of the spike. 
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Figure 1 : A representation of a sample path for the countdown process and the related membrane potential 
(case of the perfect integrate-and-fire neurons network). 



The next spike time of neuron j is modified consequently to the reception of the spike from neuron 
i and the new countdown value of neuron j after spike emission is reset according to a random variable 
denoted rjij, corresponds to the next firing time of neuron j after taking into account the reception of 
the incoming spike from i, and conditionally on the countdown value at the spike reception (see Fig 
rriPj. This random variable is called the interaction random variable. We observe that the countdown 



process is a piecewise deterministic process such as the ones described in (iDavisL 119841) . The spike 
times (ti, . . . , t n , . . .) are the only possible discontinuities of the countdown process. The countdown 
sequence is defined as the value of the iV-dimensional countdown process of the network sampled just 
after the spike emission: X n := X(t%). It is easy to show that the knowledge of the countdown process 
is equivalent to the knowledge of the countdown sequence, as we now prove. 

Proposition 1.1. There is a path-wise bi-univocal correspondence between the continuous-time count- 
down process and the discrete-time countdown sequence. 

Proof. Assume that we have a realization of the countdown process (X t )t£[o,T\- The spike times (ti, . . . ,t n , . . .) 
are given by the zeros of the process, and the value of the process just after these spike times X n := 



1 and possibly on different other variables depending on the model considered. 
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X(t+) provides the countdown sequence in a unique way. 

Conversely, assume that we dispose of the values of the countdown sequence (X n ) n£ ^. Then we 
build the related countdown process X t by "filling the gaps" using the deterministic dynamics between 
two spikes: Xq = Xq, and between two spikes, each components X t decreases linearly with slope —1. 
At time t\ equal to the minimal value of Xq, the first spike is fired, and the process X is reset to the value: 
X(t\) = X\. The same procedure is then iterated and yields the unique countdown process related to 
the countdown sequence considered. □ 

This kind of processes has already been studied in the fiel d of queueing t heory under the name of the 



hourglass model. It was first introduced b y Marie Cottrell in (ICottrelll 



1998 



and collaborators in (Flicker et al 

pool 



Cottrell and Turova, 



Turova 



19921) . stu died by Flicker, Robert 



1994) and then instantiated and studied in ( Asmus sen and Turova , 



2000 



19961) . In order to apply this modeling to the kind of 



networks of interest here, we need to define in each case the reset and the interactions random variables. 

The framework proposed here is on the equivalence of the spike times in the original network and 
the spike times emulated by the countdown process (or the countdown sequence). This property actually 
stems from the definitions of the reset and interaction random variables. We however provide a formal 
proof since this equivalence is a cornerstone of the paper. 

Theorem 1.2. Consider a network of N spiking neurons with membrane potential (and possibly addi- 
tional variables) V t . Let Vq be the initial state of the network, and assume that: 

(i) the probability distribution Yq of the N -dimensional random variable corresponding to the first 
spike times of the network starting from Vq if no interaction occurs meanwhile is known ( this 
random variable is called the initialization random variable) 

(ii) the probability distributions of the one-dimensional reset random variables (l r ' l ')j=i...jv ® n d 

(iii) the one- dimensional interaction random variables (T?jj)(i ! j)=i...jv conditionally on the last value 
of the countdown process are known and only depend on the state of the countdown process at the 
time of the spike. 

Then the probability distribution of the spike times of the original network is equal to the law of the 
spike times of the countdown process (or equivalently of the related countdown sequence) with initial 
distribution Yq, reset random variables (Y^)i=i...jv and interaction random variables (?7ij)(i,j)=l...JV- 

Remark 1. Note that the assumption that the interaction random variable only depends on the state of 
the countdown variable at the time of the spike can be weakened, and as we show below, it might depend 
on some additional variable, but needs to be independent of the underlying membrane potential. 

Proof. This theorem actually stems from the very definition of the countdown process and of the reset 
and interaction random variables. Indeed, because of assumption 1, the law of the first spike elicited by 
each neuron of the network if no interaction occurs is by definition given by Yq. Therefore, the law of 
the time t\ of the first spike elicited in the network (resp. the index of the related neuron) is equal to the 
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minimal value of Yq (resp. the coordinate index i\ where this minimum is reached), hence equal to the 
first spike time (resp. neuron index) obtained in the countdown process defined. 

The next spike time of neuron i\ if it receives no spike occurs meanwhile has by definition the law 
of Y^ ll \ which is also the probability distribution of the countdown process after the spike emission. 
All other neurons in the network will fire at a new value given by rjij by definition of these variables, 
and these are also the updated values of the countdown process, i.e. after taking into account the spike 
emission. 

Therefore, by induction, the theorem is proved: the spike times of the network are equal in law to the 
spikes emulated by the countdown process (or the countdown sequence). □ 

We note the following straightforward Markov property that follows from the specific assumptions 
of theorem [L2l 

Proposition 1.3. Under the assumptions of theorem [L2l the countdown process (resp. the countdown 
sequence) constitutes a Markov process (resp. chain). 

Note that we did not prove that the countdown process existed, nor that it was a Markov process. 
The results of theorem 11.21 and of proposition 11.31 rely on the assumption that we are able to identify 
the initialization probability distributions Yq, the reset random variables Y ( % > and the interaction random 
variables rjij. This is a more difficult task that we will deal with in sections |2] and appendices lAl and iBl in 
some specific networks where those random variables can be explicitly derived. 

One of the difficulties in event-based modeling of networks is what we call the avalanche phe- 
nomenon. This phenomenon consists in the simultaneous occurrence of infinitely many spikes. When 
connections are instantaneous and excitatory, with big enough synaptic weights, the spike emission of a 
neuron can instantaneously elicit spikes from the neurons directly connected to it, which in turn can elicit 
spikes in their neighborhood. Such a chain reaction can result in an overflow in the network, and even 
make the first neuron spike again, which closes the loop: at a given instant, the neuron fires infinitely 
many spikes. This phenomenon is an artifact of the assumed instantaneity of spike transmission and of 
the absence of refractory period in the model, this is why we will be introducing both features in our 
models in appendix|B] These features are particularly important when excitation is involved. 

2 Application to some simple types of neuron models 

In this section, we aim at building a Markov chain governing the spike times, and to this end identify the 
law of the two random variables necessary in the definition of the countdown process: the initialization, 
the reset and the interaction random variables. 

2.1 The reset and initialization random variables 

The reset random variable of the countdown process corresponds to the time to the next spike of the 
neuron that just fired a spike, if no spike occurs meanwhile. For all integrate-and-fire neurons, with 
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linear or nonlinear dynamics, when neuron i fires a spike at time t*, the membrane potential is reset at 



for the neuron to fire another spike if nothing occurs meanwhile, and therefore has the law of the first 
hitting time after t* to the threshold function 6(t) of the membrane potential process, starting from the 
reset voltage namely: 



This random variable is known in closed-form, or can be easily computed only for linear integrate- 
and-fire neurons. This is one of the reasons why we restrict our study to this class of neurons. Note 
however that, formally, this expression is valid whatever the neuron intrinsic dynamics and whatever the 
threshold chosen. The limitation only arises from the technical complexity of finding the probability 
distribution of first hitting times of stochastic processes, which remains a great endeavour in stochastic 
processes theory and still does not have any satisfactory systematic analytical or even numerical solution 
beyond linear processes. Note also that this random variable depends on t* only if the membrane potential 
process V^'(t) or if the threshold 9 are non-stationary. 

Remark 2. When the parameters of the model are not stationary, for instance when the input current or 
the threshold are deterministic non-constant functions of time, the reset random variable depends on the 
last spiking time of each neuron. Therefore in that case, in order to define the reset random variable at 
each spike time, one has to know the absolute time which can be seen as the Markov chain defined by 
the cumulative sum of the spike times. The Markov chain of the countdown process therefore involves 
the additional variable T which is updated as T n+1 = T n + inf l=1 ... N X% \ 

The initialization random variable transforms the initial condition on the membrane potential into 
the first spike time for each neuron if no spike occurs meanwhile. This random variable can be easily 
computed from the first hitting time of the membrane potential process to the threshold. Let for instance 
Vb be the initial condition of the membrane potential at time to- The initialization random variable has 
the law of the first hitting time: 



Now that these random variables are defined and their role understood, we address the case of the 
interaction random variables. 

2.2 Integrate-and-fire models without synaptic integration 

We start by considering instantaneous current-based interactions. In these cases, the reception of a spike 
results in an instantaneous current that modifies the membrane potential by an amount equal to the synap- 
tic connectivity weight, and the input current is modeled as a Brownian motion. We consider two different 
types of intrinsic dynamics: the so-called perfect integrate-and fire neuron and the leaky integrate-and-fire 
neuron. The Brownian current noise model consists in assuming that we are in the application domain of 



(i) 

a deterministic value Vy . The new countdown value at this time hence corresponds to the time it takes 





8 



the classical diffusion approximation. T his approximat i on is widely used in computational neurosci ence, 



and details can be found for instance in (Holden, 



1976 



Ricciardi and Smith, 



1977 



Tuckwell, 



1988 



2.2.1 Perfect integrate-and-flre neuron 

We start by considering the so-called perfect integrate-and-fire neuron with external input and Brownian 
noise. Between two spikes, the membrane potential of the neuron i, denoted by V^' (t) satisfies the 
equation: 



idV {i) {t) = I^(t)dt + a i dW i 



(2) 



In this equation, r, denotes the membrane potential time constant, Ie (t) the deterministic input current. 
The independent external noise sources are modeled by the N Brownian motions (WW)i<j<jv an d °~i 
are the related standard deviations. 

The equation can be integrated in closed form, conditionally to the fact that we know V® {t ) the 
membrane potential at time t = to'. 



(3) 



The neuron fires when its membrane potential reaches the threshold 9(t), after what the membrane 

(i) 

potential is reset to a constant value Vr and a spike is emitted. At this time, t , the membrane potential 
of all its neighbors j £ V(i) is instantaneously increased by the negative synaptic weight: Vw(i*) = 
VW(t-) + w i:i . 



vW(r) = e(t) 



V®(F) = V r (i) 

vU\t*) = V^(t*-)+Wi 



Assume that neuron i fires at time t* . From the result of section 12.11 the reset random variable 
of the related countdown process has the law of the first hitting time of the Brownian motion starting 
from Vr to the boundary — (8(t + t*) — f/« + * I<P(s) ds). This law can be computed using differen t 



mathematical or numerical methods (see e.g . dHolden , 



Holden 


1976 


Ricciardi and Smith, 


1977 


Plessei 


1999) 



or the review (ITouboul and Faugerasl 120071) 1. One of the most efficient methods to solve this problem 
in a general setting amounts to solving a Volterra integral equation. It is the method we use in our 
simulations when the input and/or the threshold are non-constant functions of time. Else, we will be 
usin g the well-known formulas o f the first hitting time of the Brownian motion to linear boundaries (see 
e.g. (iKaratzas and ShreveUl987h '). 

The interaction random variable is deduced from the effect of a presynaptic spike incoming at a 
synapse. If neuron j receives an inhibitory spike from neuron i at time t*, its membrane potential 
is shifted by an amount equal to Wij (see formula (01 with initial condition at time t* : V^'(t*) = 



Note that formally most of the derivations done in the present manuscript might be performed with the use of the discontinuous 
(point processes) noise Si. However, the expressions of the random variables involved in the dynamics of the countdown process 
would be harder to express in closed-form, and numerical evaluations of such random variables are generally not as efficient as the 
powerful tools of the stochastic theory. 
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V<J) + Wij y Therefore, it will never cross the threshold before time X^' (t* ), and at this time, its 
value is: 

v®(xW(f)) = e(x®(f)) +w iJ . 

Starting from this initial value the membrane potential integrates the input and the noise, and therefore 
will not spike at time t* + X^> (£*) as expected if neuron j did not receive a spike. The additional time 
after time t* + (t* ) before the next spike time of neuron j is therefore independent of the past of the 
process, and has the law of the random variable: 

m = inf {f > 0; 9{X^\t*) + t*) + Wij + J I^(s + X^(t*) + t*)ds 

+ <j l W t {l) =6{t + t* +JfW(f))}, 

i.e. the law of the first hitting time of the standard Brownian motion to the boundary 

1 f* 
in —(e(X^(t*) + t* + t)-9(X^{t*)+t*)- Wij -J I®(a + t* +X^{t*))d s y 

The interaction random variable is therefore simply T]ij(X^'(t*)) = X^\t*) + ry. The interaction 
between neurons therefore only amounts to the addition of an independent random variable on the count- 
down value. In the case of stationary inputs and threshold, this random variable has the law of the first 
hitting time of a drifted Brownian motion starting from to the constant barrier |u%|, whose density 



reads (see e.g. (|Karatzas and ShreveHl987h l: 



p^\t) = ^^ e ~ (Jf ^f L l 1R . x (t) (4) 

<Ji\2'Kt A 

In the case of stationary inputs, we can convince ourselves that the countdown process is an au- 
tonomous Markov process. Indeed, the sequence (X n ) n >o (respectively and equivalently (X t )t>o), see 
sectionQ] is a Markov chain (respectively a Markov process) with transitions given by the law of the reset 
and interaction random variables, which only depend on the parameters of the model. 

Proposition 2.1. If the input current and/or the threshold depend upon time, then the chain (X n ,t n ) n >o 
where t n is the time of the last spike fired in the network (respectively and equivalently (X t ,t) t >o) is an 
autonomous Markov chain (respectively process). 

Proof. Indeed, we have seen that in that case both the reset and the interaction random variables depend 
on the last spike time, and that the interaction random variable depends on X^'(t*). More precisely, let 
(X n , t n ) be the state of the Markov chain after n spikes (X n is the TV-dimensional countdown process 
and t n the sequence of spike times). Define i„ = Argmin i {X,l l ' ) }. The neuron i„ fires the (n + l) th 
spike at time t„+\ = t„ + Xn"\ Its countdown value X^?\ is reset by drawing from the law of the first 
hitting time of the standard Brownian motion starting from Vr to the boundary 

* ^ (V + *" + X n n) ) - I?\tn + X^ + *) ds\ , 
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and the state of the neurons j connected to i„ is updated by adding an independent random variable rji n j 
having the law of the first hitting time of the standard Brownian motion to the boundary 

t ^ — (e{X& + t n + t) - 0{X%) + t n ) - Wij - / J«(« + + t n ) ds 

°J \ Jo 



Therefore, the transition probability only involve the value of the chain at rank n, since the formulas 
only depend on t n and X n and are independent of (tk, Xk)k<n-i the previous values of the spike times 
or countdown values at previous ranks. 

We conclude from this observation and by application of theorem [L2l that the sequence (£„, X n ) n 
constitute a Markov chain the simulation of which provides access to the spike times of the network. 

□ 

In all this section we considered time dependent thresholds 0. This dependency upon time is dropped 
in the rest of the paper, since most of the classical integrate-and-fire neurons do not include this de- 
pendency. Note however that all the following developments can be performed using a time-dependent 
threshold in a same fashion as done here. 

2.2.2 Leaky integrate-and-fire models 

The model We now take into account the leak of the membrane potential but no synaptic integration. 
The membrane potential is governed by equation: 

ndV^ (t) = (-V« (t) + (t))dt + <TidW t (i) 
V«(*~) =0=*.y«(t) = v r w 

where the (W7)i<t<jv are independent Brownian motions. Therefore, the membrane potential yMM 
in the absence of correlation is an Ornstein-Uhlenbeck process. Between two spikes, the membrane 
potential can be written as: 

ft rt 



V®(t) = V®(t )e- < -*- to)/T *+ e {s ~ t)/ ^4 l \s)ds + <i l e (s - t)/T * dW® (5) 

Jto Jto 

together with the spiking condition: 

V«(*) =Vr {i) 

V® (t) = (*") + Wij Vj £ V(i) 



VW(t~) = 



Let t* be the time when the neuron j receives a spike from neuron i. We define := V^(t*~) and 
X ( J ] := XW(f-). We have: 

V {j) (t* +t) = (Vy } + loyje - */^ + - / e (s -*) /T '/W(s + t*) ds 



J JO 



+ (* ei'-tVn aw® (6) 

T 3 JO 
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Let the membrane potential of the neuron j if it did not receive the incoming spike at time t*. It 
satisfies the equation: 

VW(t* +t) = V^V'/^ + - f e(-*V r >lW(s + t*)ds + ^ f e ( s ~*)/^ dW® 

Tj Jo Tj Jr. 

Therefore we have: 

V u) (t* +t) = V u} (t* +t) + u^e - */ 7 * 

(i) 

(note that this property simply stems from the linearity of the LIF neuron). For t = XI we have 
VW (t* + = 9 and, from ©, we have: 

vWp+xP+t) = (9+w tJ e- x * ) ^)e- t ^+- f e^/^ 1^ (s+t*+xi j) ) ds+^ f V*"')/^ dW W 



We therefore deduce that the next spike time, i.e. the next hitting time of the barrier 6 by the process 
conditionally to the random variable X* is the sum of X* and an independent random variable 
whose law is equal to the hitting time to the barrier of the process (0 with initial condition V^' (0) = 

« + iD,je"^ 3)/T = and with the time shifted input current I { e 3 \t) :=I^\t + t* + xi 3) ). 



Tjij - inf {t > 0; U u) (t) = 9\U U) {0) = 9 + w^e-**' /t > } 



where U^'(t) is the solution of equation <j3j with the specified time-shifted current. 

The problem of the first hitti ng time of the LIF neuron with constant or curved boundaries was ad- 



dressed in many articles (see e.g. (lPlesserlll999tlTouboul and Faugeras 



2007) for a review and references 



therein), no closed-form solution is available, and this distribution can be efficiently computed by solving 
a Volterra equation, for instance. 

An important remark is that this random variable only depends on X^\ Conditionally to xi^ , the 
added random variable is independent of the past of the process, so the sequence X^ is Markovian. 
Furthermore, the network countdown process dynamics is autonomous: we do not need to refer to the 
underlying membrane potential to describe its evolution. This is very interesting since we can study and 
simulate this random variable by itself. Therefore, the variable (X t ), possibly augmented with the time t 
if the input current or the threshold are not stationary, is a Markov process, and this process sampled at 
the times of the spikes is a Markov chain. Furthermore, the law of the zeros of this process is equal to 
the one of the spikes of the underlying network. 



Event-based simulations The Markov chain we built provides a very natural event-based simulation 
algorithm for inhibitory stochastic networks. The method consists in simulating the Markov chain de- 
scribing the time of the spikes for each neuron. We have seen that simulating the times of the spikes was 
equivalent in law to simulating the membrane potential, from the spikes viewpoint (see theorem [L2l . The 
event-based simulation consists in building this Markov chain. Simulating this Markov chain requires to 
draw at each spike time from the laws of the reset and the interaction random variables. We have seen that 
these random variables can be expressed in most of the cases as first hitting times of random processes. 
In the cases where these laws are known in closed form, a very efficient simulation procedure can be 
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used. If it is not the case, then we will have to evaluate these random variables. Assume that at the initial 
time to the values of membrane potential of each neuron and of the additional variables of the model are 
known. The initial countdown value for a given neuron will be simply computed as the first hitting time 
of its membrane potential process to the threshold, starting from this initial condition and can therefore 
be computed in the same way as the initialization variable. From this initial time, the principle of the 
algorithm is to build the discrete-time Markov chain containing as a variable the countdown process that 
gives the times of the spikes (we have seen that sometimes additional variables were necessary). Then 
to deduce the state of the chain at time n + 1 knowing its state at time n, we use the recursion relation 
described in previous sections. 

• We first identify the neuron having the lowest countdown value, which amounts to finding the 
minimal value in a list of N elements, an elementary operation efficiently coded. This neuron is 
the one that elicits the first spike. 

• When this neuron is identified, we move the clock forward to this time, and draw the new state 
of the network: the neuron that just fired a spike is reset by drawing from the law of the related 
reset variable and the other neurons states are updated by drawing from the law of their respective 
interaction variables. Once the states of all neurons have been updated, the simulation proceeds. 

Let us for instance compare the natural Monte-Carlo simulation and the event-based simulation, in 
terms of precision and computational efficiency, in the case of a simple two-neurons perfect integrate-and- 
fire network model with instantaneous currents, each neuron having a constant threshold 6i and receiving 
a constant input I{ . In that case, based on the previous analysis, the reset and the interaction random 
variables involved in the dynamics of the countdown chain respectively have the following densities: 




-v* / {Qi-vl-ii tf 



1^1 exp ( - {w ^ +h t] ' 



(7) 



In figure|2]we compare the results of the simulations of the two algorithms on 500 000 realizations of the 
network and observe as predicted a very good match between the two distributions. In order to directly 
simulate the sample paths and the spike times, Monte-Carlo methods prove relatively inefficient, since 
the error in th e spike time is of the order s/dt if dt is the time step of construc tion of the sam ple path (see 



(IGobetl 120001) '). In order to improve the accuracy of the algorithm, Gobet in (IGobetL 120001) proposed the 



interesting idea to interpolate the threshold crossing using the law of the maximum of a Brownian bridge 
conditioned on the values it takes at the discretized times. For the event-based algorithm, we tabulated the 
functions given in equations © for generality^. We use both natural and Gobet's refined Monte-Carlo 
method and compare the precision and the execution time of the numerical algorithms. Implemented in 
Matlab®R2008b, and run on a MacBook 2Ghz Interl Core 2 Duo with 2GB 1067 MHz DDR3 memory, 
Gobet's enhanced Monte-Carlo algorithm for 50 000 sample paths runs in 716.68 seconds. On the same 



J The case of PIF neurons is the only case where both the interaction and the reset random variable can be expressed in closed 
form with usual functions. 
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(a) Spike time distributions, neuron 1 



(b) Spike time distributions, neuron 2 




1 2 „ 3 
time 




(c) Spike time distributions, neuron 1 



(d) Spike time distributions, neuron 2 



Figure 2: Spiking probability as a function of time for a fixed initial condition: histogram computed over 
500 000 realizations of the process. Com parison of the spike time distributions simulated by a Monte- 



Carlo method (green solid line), Gobet's (IGobetL 120001) Monte-Carlo improved method (blue solid line) 



and by our event-based simulation method (dashed red line). Left: Neuron 1, Right: Neuron 2, top row: 
symmetrical network with input current h = h = 1, reset V} = V? = 0, thresholds 9% = 62 = 1 
and connectivity weights W12 = W21 = —0.2. Bottom Row: Asymmetrical network, input current 
I\ = I2 = 1, reset — = 0, thresholds 6\ = 1, 62 = 1-3 and connectivity weights W12 = —0.1 
and W21 = —0.5. Note that because of the prohibitive execution time of Monte-Carlo methods, these are 
coded in C++ for the purpose of this simulation over a large number of sample paths. The time interval 
chosen is [0, 4] and the time step for Monte-Carlo algorithms is dt = 0.01. 
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machine, we simulated the natural Monte-Carlo algorithm using the same discretization time-step, in 
which case the precision is not at the level of both Gobet's and the event based algorithm as we observe 
in Figure |2] the execution time is of 19.04 seconds. Note that in order to achieve the same precision 
on the spike times as Gobet's method with discretization time step dt, one should use a time step of the 
order of dt 2 for the regular Monte-Carlo method, in which case the execution time of the event-based 
algorithm for the same parameters shall be a hundredfold larger, a prohibitive execution time for such 
a simple model in such small dimensions. On the same machine and with the same number of sample 
paths, the event-based algorithm runs in 2.46 seconds. Simulating the event-based algorithm in this case 
therefore yields a gain of the order of 30. It therefore clearly appears that the event-based algorithm is 
faster. Moreover, increasing the precision of the latter only involves a more precise tabulation and does 
not affect the efficiency of the algorithm, whereas dividing for instance by 10 the time step of a Monte- 
Carlo algorithm multiplies by 10 the execution time. Note also that the number of neurons in the network 
and the time interval over which the simulation is performed will not modify this ratio. 

Let us now describe the result of figure Fig. |2] In the top row we show the case of a symmetrical 
network. In the absence of noise and interaction, each neuron fires with period 1 . The presence of noise 
results in the fact that the distribution of the first spike time will have its median value at 1 with a certain 
spread around this value. Therefore one of the neurons (say neuron 1 without loss of generality) will fire 
first. It instantaneously inhibits neuron 2 that will now fire after some delay. Because of the symmetry of 
the network, this phenomenon accounts for the multimodal nature of the distribution. Consecutive spikes 
reproduce the same scheme. The asymmetrical case can be understood in the same fashion: neuron 1 
fires with a median ISI equal to 1, and neuron 2 with median ISI 1.3. Neuron 1 fires more frequently, but 
its connectivity weight is smaller and affects less the state of neuron 2 than neuron 2 affects the state of 
neuron 1 . These phenomena result in shifting the peaks of emitted spike times accordingly. 

Therefore, we have verified that in the case of instantaneous synapses, the sequence of spike times 
constitutes a Markov chain, whose transition probability was explicitly derived from the law of the first 
hitting time of the membrane potential process to some boundary. This result extends to cases where 
on takes into account synaptic integration as we now show. In appendixes [A] and [B] we show that these 
results can be extended to more complex cases, but fail to do so for conductance-based models. 

2.3 Integrate-and-fire models with exponentially decaying synaptic currents 

In this section we consider that the spike integration and the noisy integration are not instantaneous. Note 
that by exponentially decaying synaptic currents, we are only concerned with the noise-generated current, 
and we will consider in this section both instantaneous and exponentially decaying spike integration. We 
start by considering, similarly to the previous case, the perfect integrator case first before addressing the 
case of the leaky integrator. 
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2.3.1 Perfect integrate-and-flre model 



A different way for improving the perfect integrate-and-fire model is to consider that the inputs received 
by the neuron i are integrated at the level of the synapse with a characteristic time constant ri 1 ' 7^ 0. In 
that case, the membrane potential satisfies the equations: 



whose solution reads: 



/W(u)du + r 8 w (l 



J«(i) + /«(i) 
-/ S W (t)dt + aidW^ 



e _( w _„)/ T w dT ^i) dv ^ (8) 



to ^to 



or in the leaky integrate-and-fire case 

(t) = - V® {t))dt + 7 e (i) (i)tft + rf } 1 < i < N 



rPdl®® = -I*f ) {t)dt- 
in which case the membrane potential reads 



(9) 



GidWl 



V il) (0) + / 1^ (s) ds + rj 4 ' (1 - e-'/ r = W ) J« (0) + , 



-(s-u)/t^ dw (i) d 



and a spike is emitted when the membrane potential reaches a threshold 9 which we assume constant. 

We now consider the related countdown process (Xt)t and compute the random variables that are 
necessary to the definition of its dynamics. The related reset random variable in the case of the expo- 
nentially decaying synaptic currents has the law of the first hitting time of the Doubly Integrated Process 
(DIP) J* r o " e -(-"-«)/^ w dwi l) dv or Integrated Wiener Process (IWP) /„ du to a curved bound- 



times canbe evaluated based on methods proposed e.g. in (iMcKean , 



(McKean 


1963; 


Goldman, 


1971; 


Lachal 



1991 



1996b i n the case of constant input a nd constant threshold, or approximated using the framework 



developed in (ITouboul and Faugeras , 



2008). Note that in this latter case no closed-form solution can be 
provided except in the case of the IWP with a piecewise quadratic polynomial input. 

Let us now compute the interaction random variable, in the cases of instantaneous or exponentially 
decaying synaptic spike-integration currents. 



Instantaneous synapses In the case of instantaneous synaptic integration of the spikes, the linearity of 
the equation yields the fact that the interaction random variable has the law of the first hitting time of the 
threshold 9 of the membrane potential process starting from (9 + ;e -X * J, / T i 3) , li l \xi^))) to reach 
the threshold 9, and can therefore be computed using the same approximations or formulas depending 
on the model we choose and the type of input current considered. Furthermore, because of the strong 
Markov property of the membrane potential process, the interaction random variable is independent of 
the trajectory of the voltage and only depends on Xj?' . 
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Synaptic integration Taking into account synaptic integration in the connections makes more sense: 
it corresponds to the assumption that whatever its origin, the inputs received at the level of the synapse 
are integrated in the same fashion. In that case the effect of an incoming spike on a postsynaptic neuron 
is added instantaneously to the synaptic current. Therefore, using the same technique as before, we can 
obtain the law of the interaction variable. For the perfect integrate and fire neuron with exponentially 
decaying synaptic current, if neuron j receives an incoming spike from neuron i at time t*, the law of 
this random variable is deduced from the law of the first hitting time of the related DIP starting from 
(d + w tj T { s 3) (1 - e- x * 3) )/ T "' ), I s (t* + Xi 3) ^j + w lj e- x * ) \ to reach the threshold 9. In the case 
of the perfect integrate-and-fire neuron with perfect synapses, the law of the interaction random variable 
is given by the law of the first hitting time of the related IWP starting from (6 + xi^Wij , I s (A* ) + my ) 
to reach the threshold 9, and this random variable is independent of the membrane potential and only 
depends on a* . 



Building a Markov chain In these two cases the countdown process is not a Markov process. If 
we consider in addition the value of the synaptic current I s (t n ) at the time t n , we obtain a Markov 
chain modeling the times of the spikes. More precisely, consider (X n , I s (t n ),t n ) the countdown process 
together with the synaptic current at the next spike time and t n the last spike time. The first spike will 
be fired by the neuron i n having the lowest countdown value. It will fire at time t n +i = t n + Xn^ ■ 
Its countdown value and the value of the synaptic current at the next spike time is computed by drawing 
from the law of an independent p air composed of the next spik e time and the relative location of the 
synaptic currents at this time (see (ITouboul and FaugerasL 120081) ). The countdown value and the future 
synaptic current of each neuron j E V(i) is updated by drawing from the law of the first hitting time of 
the membrane potential of neuron j, and the other neurons states are unchanged. This chain is a Markov 
chain whose spike times are equal in law to spike times of the original network, by direct application of 
theorem [L2l and of proposition ! 1.31 



2.3.2 Leaky integrate-and-fire model 

Our last example features both a leaky intrinsic dynamics and synaptic integration of the noise. In this 
case the membrane potential and the synaptic noise are coupled via the differential equation (0 with the 
spiking condition given by: 

V^{t-) = 9^{ ... W r (10) 

Qualitatively, when a spike is received by a neuron, the synaptic current I s integrates the spike and 
the effect on the membrane potential is smoother. Therefore in this model it is interesting to consider 
post-synaptic pulses having the same dynamics as the noise integration, i.e. solution of the differential 
equation: 

= _ I (i) (t ) 
s dt a [ ) - 
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The very same analysis could be done if we considered an instantaneous spike integration. We omit it here 
since it follows in a straightforward fashion from the upcoming analysis of a slightly more complicated 
case. 

The reset random variable is given by the first hitting time of the membrane potential stochastic 
process. The same type of calculations as in the previous section yields the relation: 

V u) (t* +t) = V u) (t* +t) +e~ t ' T iw i j— (11) 

a>j 

for j e V(i) and Tj ^ ri . ctj = -^y — — and again V^\t* + t) is the membrane potential of 

s (7) 
the neuron j without any interaction. We can see that after the time XT , the membrane potential of j 

is 9 + Wije^ t / Ti 1 ~ e - 3 . The evolution of the potential V^' after t* + X* and conditionally on 

X* and I s (t*) is independent of the past, so we have to wait for the process dTOb to reach the threshold 

from the initial condition 6 + Wije~ t ' Tj 1 ~ e - 1 and with the "time and space" shifted currents 



-4 J) 

li(t) := I((t + t* + X* W ) + w^e'^P 7 . 

In the case Tj = rj we replace the expression 1-6 a 3 by X\^ Wije~ l l Ti , and the change in the 
currents is the same. 

Therefore, the variable (X t ,I s (t)) is Markovian and we deduce the precise firing times from its 
study. This Markovian variable requires to evaluate the law of the first hitting time of a DIP to a curved 



boundary, which can be achieved using the technique developed in (ITouboul and Faugeras , 



2008). As in 



the case of the perfect integrate-and-fire neuron, we can show that this process, possibly augmented with 
the time process t, and possibly sampled at the times of the spikes, satisfies the Markov property and the 
law of the zeros of the Markov process is the same as the law of the spikes of the underlying network. 

In the case of non-stationary inputs, the dynamics of the countdown process together with the synaptic 
current at the next firing time and the last spike time can be described as follows: let us note (X n , I s ,ni t n ) 
the previous Markov process after the n th spike. The next spike will be fired by the neuron i n having the 
lowest countdown value. It will fire at time t n+ i = t n + X„ \ Its countdown value will be reset to the 
first hitting time of the related DIP, which can be evaluated as before. Therefore, by drawing from the 
law of this pair, we have the new countdown value and the future synaptic current at the time of the next 
spike for i n . Similarly, each neuron j E V{i n ) is updated according to the law of the first hitting time 
of the related membrane potential starting with an input current given by Is% to reach a given threshold, 
and therefore the new countdown value and the future location of the input current are computed at the 
same time using the same results as before. The other neurons states are unchanged. It is clear that the 
law of the spikes is the same as the law of the zeros of the countdown process (see theorem [L2b . 



3 Excitatory Networks 

We have shown in the previous section that when the network is inhibitory, the strong Markov property 
of the membrane potential process and the linearity of the equation allowed to define an event-based 
description of the network activity. In appendices [A] and |Bj we show that these results can be extended 
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to other more complex cases, but not to conductance-based models. In this section, we are interested in 
extending these results to the case of excitatory interactions. 

As already mentioned, a well-posedness problem arises when taking into account excitatory interac- 
tions. This is due to the instantaneity of the interactions which might lead to avalanche phenomena, i.e. 
self-excitation of the cell through the network at a given instant. A suitable way to avoid such problems 
is to take into account two important biological phenomena: the synaptic transmission delays and the 
refractory period. Indeed, both phenomena prevent the occurence of avalanches by upperbounding the 
firing frequency and forbidding instantaneous self-excitation. We show in appendix [B] that they can be 
introduced in the model and only amount to transforming the interaction random variables involved in 
the simple case where no delay or refractory period is taken into account. This is why, in the case of exci- 
tatory connections, we focus on the arrival of a single incoming excitatory spike at a given time, keeping 
in mind that, thanks to a suitable implementation, these two phenomena can be taken into account. 



3.1 General facts 

We start by formally presenting the method on an unspecified general linear model of the same type as 
the ones covered in section [2] and some of the models of appendix lAl where the Markovian description 
of the countdown process is possible. Let us assume that at time t = t* , neuron j receives an excitation 
coming from neuron i. Note that, strictly speaking, this scenario covers cases including synaptic delays 
and refractory periocQ. 

The general principle of our approach for defining the future spike time of neuron j is based on 
characterizing the value of the membrane potential (and of the additional variables needed to define the 
state of the membrane potential involved in its dynamics) of neuron j at time t* . Let us first treat the 
unidimensional cases of the PIF and LIF neurons with instantaneous synaptic currents. Let us assume 
for that both the last spike time ii ast and the value of the membrane potential (ti as t) at this time 
are known. Then the law of the countdown value after spike reception can be written, using the total 
probability formula: 



P 



xv\e+) g dt\x^(t*),t last ,v^(t Usi ) 



p 



j T 8 t ., Vj+Wi . et* + dt\v^(t*) = iJph*j)(f) G d Vj \X^(t*),t last ,V^(t lm ) 



(12) 



where ■ 7 rf u is the first hitting time of the threshold of starting from initial condition u at time s. 
We used the traditional notation in the stochastic literature P(£ g dz), which represents the probability 
measure of a random variable ( at z. This can be seen as the infinitesimal probability of ( to be in dz. In 
all the cases we consider the random variables have densities with respect to the Lebesgue measure, and 



^Indeed, if the action of the spike fired by neuron % at time t; affects neuron j after a delay Ajj, then t* = if + Ajj and if 
tj is the last firing time of neuron j and the refractory period corresponds to modulating the interaction by a function k, then the 
synaptic weight Wij is multiplied by n(t* — tj). The way to include such phenomena in a generic manner is presented in details 
in appendix FBI 
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this quantity is equal to P(;(z) dz where is the probability density of £ with respect to the Lebesgue's 
measure dz. 

Formula (fT2l readily extends to multidimensional cases such as the bidimensional PIF and LIF neu- 
rons with synaptic integration treated in this manuscript. Indeed, assuming again that the last firing time 
tiast and the values of the membrane potential and of the synaptic current at this time (V^ (ii as t), Is (t\m)) 
are known, the total probability formula yields for the probability of the next firing time: 



P 



X«(t*+) g ^|X«(t*),t last ,(y«(t last ),7^)(t last )) 



P 



x P 



V^(t*) G do s ,lP(t*) G ^|X«(f),t last ,y«(t last ),/«(t last ) 



(13) 



The first term in the integral of equations (fT2l ) and ( fT3l involves the probability distribution of the first 
passage time of the membrane potential conditionally to a fixed initial value, which corresponds exactly 
to the problem of finding the initialization random variable. As shown in section |2~T1 this variable is 
directly related to first hitting times of the membrane potential process, and different ways for computing 
this law are discussed depending on the model considered throughout section [2] The second term 
involves the probability distribution of the membrane potential (and possibly additional variables such as 
the synaptic current) conditioned on an initial state and on the value of its first hitting time to the threshold 
boundary. This law can be computed through the use of Bayes' formula for conditional probability. Let 
us treat for instance the case of a one-dimensional model (the general case is treated in a totally identical 
manner, but involves more intricate formulae). We denote for the sake of compactness of notations r the 



first hitting time 3 t 



t*,v«(t«) 



of the process (Yt )t>t* to the threshold 9, with initial condition at time 
t* equal to V^(t*). We have, using Bayes' formula for probability densities and the Markov property 



P 
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= t* + X ( t P 







(14) 



The formula (fT2l , thanks to the relationship (fl4l only involves quantities that are already computed 
or easily derived. Indeed, it involves i) the law of the first hitting time of the voltage process to the 
boundary conditioned on its initial datum, problem which is exactly solved as the computation of the 
reset random variable in section 12.11 and ii) the probability density of the Markovian transition of the 
voltage potential between t\ m and t* . This transition is very easy to estimate from the properties of the 
underlying voltage process. Indeed, in all the cases treated, the voltage potential is a Gauss-Markov 
process, hence the transition probability is a Gaussian process whose mean and covariance matrix are 
very simple to compute as a function of the parameters of the membrane potential process. This implies 
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that the characterization of the law of the first hitting time of the membrane potential to the spiking 
threshold allows characterization of both the inhibitory and the excitatory interaction random variable. 

Formulae (fT2l and (Tl3T > are compact formulations for the law of the next spike time of neuron j. It 
has an intricate form involving an integral that might be difficult to compute. However, we can express 
the interaction random variable for excitatory interactions in a more intuitive and efficient way using to 
use a two-step algorithm, providing a straightforward simulation method. Indeed, in order to compute 
the nest spike time, one only needs to draw a value of the next firing time from the law given by dTZb or 
( Tol l. And to complete this step, the computation of the integral involved in the definition is not necessary, 
since drawing a realization of this random variable can be more easily performed as follows: 

• first, draw the value Vj (or (vj, ij)) from the law of the membrane potential process (and possibly 
additional variables) conditioned on initial datum and on the value of its first hitting time, com- 
puted through Bayes' formula on conditional probabilities ([Pil l, or its multidimensional versions 
depending on the neuron model, 

• second, draw the value of the first hitting time of the diffusion process V' J to the boundary starting 
from initial condition Vj + Wij (or (vj, ij + Wij)). 

If the probability distribution of the first hitting time of the membrane potential to the spiking thresh- 
old is well characterized, both random variable involved is also well defined. The characterization of such 
first passage time problems was treated in section [2] and therefore allow event-based description of the 
spike times. We now make explicit the transition density in the cases treated in section |2] and emphasize 
on the case of the PIF neuron since its simplicity allows compact and clear presentation. Other cases 
are less detailed since the analysis of the system and the Markov property are derived in a very similar 
fashion. 

3.2 PIF Neuron with Instantaneous Interactions 

Let us start by dealing with the case of the PIF neuron with instantaneous interactions. The dynamics 

of the membrane potential is given by equation (fJJ, in which case each membrane potential is an inde- 

(i) 

pendent Brownian motion. Let to denote the last spike time fired in the network and by V the value 
of the membrane potential of j at this time. Neuron j is expected to reach the threshold 9 and fire at 
time ti = to + (to)- Meanwhile, at time t* G [to, ti], the neuron receives an incoming spike from 
neuron i that instantaneously shifts the membrane potential value by an amount equal toy > 0. Bayes's 
formula (TPfl i allows characterizing the value of V^\t*) of the membrane potential of neuron j at the 
time it receives the excitatory spike from neuron i. We have the following: 

Proposition 3.1. A network of perfect integrate-and-fire neurons with excitatory and inhibitory connec- 
tions has a Markovian event-based description. 

Proof. Let us denote as usual by (t n ) n ^ the sequence of spike times fired in the network and by 
V n = (V^ ,k = 1 . . . N) the values of the membrane potential of all the neurons at each time. We 
consider the sequence Z n = (t n ,V n , X n ) n where X n is as usual the A^-dimensional countdown process. 
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Let us assume that the state of the chain Z n is known. We aim at deriving the transition to the 
following state Z n +i. Let us denote by i n the index of the smallest coordinate of X„. It is the index 

(i ) 

of the neuron that will fire first. This neuron will fire at time t n+ i = t n + X„ and its countdown 
value X^?\ can be computed by drawing in the law of the reset random variable. is simply equal 

to V^ n) . Let us now consider the effect of the spike fired by i n on one of its neighbors, say j. If the 
interaction is inhibitory, we can apply the results of section|2]to compute the interaction random variable. 
If the interaction is excitatory, thanks to formula (fT2l , we know that provided that we are able to define 
Ki+i' we can derive the interaction random variable. 

( i) ( i) 

The law of V^li is uniquely determined by the value of Vn at time t n and the value of the first 
passage time after t n , namely t n +Xn , because of formula (TPfl i. and therefore only depend on quantities 

(i) 

given by the state Z n . The knowledge of V^V 1 eventually allows to compute the probability density of the 
next spike time of j by drawing in the law of the first hitting time of the underlying membrane potential 
process to the spiking threshold, allowing eventually to compute X^__j_ 1 . 

From our knowledge of the multidimensional variable Z n , we have exhibited a way to draw the next 
state Z n+ i of the chain, regardless of the values of Zk for k = 0, . . . , n — 1, implying that the sequence 
(Zk)k enjoys the Markov property. 

□ 



3.3 LIF with instantaneous synapses 

In the case of the leaky integrate-and-fire neuron with instantaneous synaptic integration, a similar anal- 
ysis leads to the fact that the new spike time after interaction has the law of the first hitting time of the 
membrane potential process to reach the spiking threshold 6 conditionally on the fact that it has the value 
V^\t*) + Wij at time t*. Here again, Bayes' rule allows computing in a simple form the distribution 
value of VV)(F) as a function of the law of the first hitting time of the membrane potential process. 
Based on this law, it is very simple to extend the proof of proposition |3.1| to the present case and therefore 
provide a Markovian description of the spike times for this kind of network models. 

3.4 Exponentially decaying synaptic interactions 

In the case of exponentially decaying synaptic conductances, the same framework allows deriving the 
interaction random variable in the excitatory case. Let us consider the case of the perfect integrate-and- 
fire neuron with synaptic integration and excitatory synapses. In that case, the interaction random variable 
is given by equation (fT~3T > and involves the value of the membrane potential at the time the postsynaptic 
neuron receives a spike, which is evaluated using Bayes' formula in the same fashion as we did for the 
PIF neuron in <fT~4~b . 

More precisely, assume that the neuron j receives an excitatory spike from neuron i at time t*. 
Provided that this value of the membrane potential V^(t*) and of the synaptic current I^\t*) are 
known, the next spike time fired by neuron j has the law of the the first hitting time of the spiking 
threshold 9 of the membrane potential process with initial condition (Vj(t*), Ij(t*) + Wij) at time t*. 
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The law of the (V^(t*),Ig (t*)) is derived from Bayes' rule in a similar way as we did for the PIF 
neuron (formula (O). This allows defining a Markov chain containing the times of the spikes, in a very 
similar fashion as done in proposition l3.ll 

The case of the leaky integrate-and-fire neuron with exponentially decaying synaptic currents can be 
analyzed in the same manner. 

Conclusion 

We have developed an event-based mathematical framework for the study of stochastic integrate-and-fire 
neural networks. We show that in the case of the linear models treated, the spike times can be defined 
through a Markov chain, a method which provides an elegant way to describe the spike times and an 
efficient way of simulating the network. 

The approach is quite versatile, since it can deal with a large class of neuron models and different 
types of interactions, i.e., it can cope with transmission delays, absolute and relative refractory periods. 
However, it is necessary to assume that the neuron membrane has a linear free dynamics in order to 
be able to find the transition probabilites in closed form, as a function of the first passage-time of 
the underlying membrane potential process to the spking threshold, which is related to the law of the 
interspike interval of an isolated neuron. These variables as always well defined, in most of the cases 
treated, either closed form formulae or numerical methods are available, allowing efficient event-based 
simulation. This is for instance the case of the PIF neuron with constant current, where a closed-form 
formula for the first passage time to a constant boundary is known. 

One of the strengths of the approach is that is makes use of the probabilistic nature of the firings. 
More precisely, attempts to map the actual spike times of the voltage process to a realization of a Markov 
chain governing the times of the spikes are out of reach. As noted in the excitatory case where an explicit 
dependency on the value of the membrane potential process appears, the actual spike time for a given 
realization of the voltage depends on the specific sample path produced. However, by releasing the 
constraint of pathwise agreement, we proposed a well-defined framework based on the notion of equality 
in law which is the correct notion of equality in a probabilistic framework. In other words, instead of 
mapping all the sample paths of the voltage values to one sample path of the countdown process, our 
construction builds a well-defined Markov chain which exactly reproduces the distribution of the spike 
times of the original network model. The prior to this approach is highly relevant to the underlying 
biological phenomenon: it consists in considering that it is the probability distribution of the spike times 
that contains the information of the network. 

We were able to produce expressions for the transition matrix of the spike time variables in some 
specific cases, something which had not been done before, and were able to address most of the cases 
that were addressed in the deterministic case. However, all our models remain somewhat over-simplified. 
The extension to a wider class of neuron models including more realistic features, in particular nonlin- 
earities and conductance-based interactions would be highly desirable. As shown in appendix |Aj the 
extension to conductance-based noise or interactions is a very complex task. The extension to nonlinear 
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neuron models is also difficult. In determ inistic networks, the only exception is the work of Tonnelier 



and collaborators dTonnelier et al 



2007) with quadratic neurons making use of the invertibility of the 



solutions, known in closed form. In the stochastic case, there is no closed-form solution to the membrane 
potential (now stochastic) equation which prevents from using similar methods. But an even deeper and 
more fundamental problem is the fact that the reset random variable will involve first hitting times of a 
quadratic process of type: 

dX t = Xfdt + adWt. 

The problem of characterizing such hitting times is still an active area of research in mathematics and no 
solution is currently in view. Moreover, as we mentioned above, the knowledge of the first spike time 
if nothing occurs meanwhile still implies that one has to cope with a complex correlation structure that 
makes the law of the interaction random variable quite complicated. 

The framework we have developed is only concerned with exact simulation of the spike times in 
a stochastic network. The exact simulation has been proved possible in few simple cases, mostly the 
same as in the deterministic case. However, we observed that our take on this problem did not allow 
generalization to more complex neuron models. This observation suggests to relax the constraint of exact 
simulation and develop approximate simulation frameworks based on spike times. The development of 
approximate event-based description with a controlled error on the spike times is highly desirable and 
the difficulties encountered in our study when dealing with exact simulation suggest further exploration 
along these lines. One of the possible ways to approach this problem would be to approximate the initial 
dynamics by one of the models that can be treated exactly, and for this purpose to find regimes and/or 
scales in which neurons can be approximated by linear integrate-and-fire models. The main difficulty is 
to deal with the complex correlation structure of the different random variables involved in this putative 
event-based framework which can be the source of an accumulation of errors. 

The Markovian model developed in this article gives a compact and simple way for describing the 
complex dynamics spiking networks. Besides these advantages, it is close to some renewal processes 
with interaction that have been studied in queueing theory. Though the results obtained in this domain 
do not apply so far, extensions of these theories might be reachable to prove properties such as the 
ergodicity of the models. Indeed, the theory developed only accounts for constant interaction random 
variable, i.e. whose law does not depend on the st ate of the networ k, and that have a finite expectation. 



This is also a very active field of research (see e.g. dTouboul , 



2008bl) ). Our opinion is that this framework 



allows to use the powerful tools of communication networks theory to analyze the properties of several 
important models of neuronal networks in a potentially mathematically more tractable way than the study 
of the stochastic membrane potential. This would also open the door to the mathematical study of the 
macroscopic behavior of large networks using the hydrodynamics limits d eveloped for large qu euing 



processes, to infer and model collective behaviors of such networks (see e.g. (iGromoll et al 



2008)). The 



simulation efficiency of the algorithm described he re opens the way to very large scale simulations of 
networks, a current active area of research (see e.g. (Ilzhikevich and Edelmanl 120081) ), and would allow a 
numerical study of finite-size effects in this type of networks. 
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A Inhibitory Networks for a wider class of integrate-and-fire mod- 
els 

We study more refined descriptions of the neuronal activity and show that in these cases an event-based 
description of the network activity may fail. 

A.l LIF model with general post-synaptic current pulse 

We consider a LIF neuron described by (0, but consider that the spike integration is governed by a 
general postsynaptic potential function a( ), the Green function of the differential operator C. In that 
case, the synaptic input to the neuron j, denoted Ia\ is solution of Cla = between two spikes, and 
when neuron i fires a spike, the synaptic current I a is added the synaptic weight Wjj . In the case of the 
LIF neuron, the equations of the model read: 

Tl dV t {l) =(-V t + I e (t) + ti ] (t)) dt + a dW® 
CljP(t) =0 

When the membrane potential V^ 1 ' reaches the threshold 9, a spike is fired, the membrane potential is 
reset to the value Vv and the synaptic current I a is added the synaptic weight Wij . The reset random 
variable is therefore simply the first hitting time of the Ornstein-Uhlenbeck process to the threshold, with 
synaptic current I a (t) = 0. If neuron j receives a spike from neuron i at time t*, then the synaptic current 
I a is added . Therefore, the updated membrane potential at time t* + X (t* ) reads: 

9{t* +I w (f)) + ™, je - x(,l ( f ) / a(s)e s ^ds (15) 

Jo 

and Ia\t* + X^lt*)) is the value at time t* + X^(t*) of the solution of the ordinary differential 
equation Lii = with the initial condition at time t* equal to Ia (t*~) + Wy . Note that if the operator 
C is of order d greater than one, the initial condition involve the derivatives of I a ^ up to order d— 1. These 
are unchanged by the spike reception. Hence the additional time induced by the reception of a presynaptic 
spike from neuron i has the law of the first hitting time to the boundary 8 of the stochastic process 
with initial condition at time t* + X^(t*) given by O and I a j \t* + lW(f)) + w tj W (f* )) 
where ^ denotes the flow of the linear equation = 0. This random variable has a law independent 
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of the value of the membrane potential at time t* and therefore one can build a Markov chain governing 
the times of the spikes. In details, in order to take into account this synaptic integration of spikes in 
our framework, we have to extend the phase space of our Markov chain. The Markovian variable we 
consider is the process (X t , I a (t))t>o- When a neuron i elicits a spike, i.e. when its countdown reaches 
at time t* , its countdown value is reset by drawing from the law of the first hitting time of the membrane 
potential with initial condition (V r ,Ia(t*)) to the threshold and for all neuron j € V(i), their spike- 
induced current Ia\t*) are instantaneously updated by adding the synaptic efficiency Wij : Ia\t*) = 
Ia^ (t*~ ) + Wij . Simulating this Markov process, that can be sampled at the times of the spike emission, 
is equivalent from the spikes point of view as simulating the whole membrane potential process (see 
theorem 1 1.2b . 



A.2 LIF models with noisy conductances 

We consider the case of noisy conductance-based models. In this case, the membrane potential is solution 
of equation: 

dV® = {I^\i)-\(V t {t) -V r { P v ))dt + I^\i)dt + a l g l {V t {r> -V^dWl 

r-\ ' (16) 

VW(t-) = 0=> V«(t) = V r {l) 

(i) 

The term Is corresponds to the current generated by the reception of spikes emitted from neurons in the 

H) 

network. V re v is the reversal potential of the synapse. 

Between two spikes, the stochastic differential equation can be integrated in closed-for m. For in- 



fo 

stance, provided that V to is known, we have for t > t a (see e.g. (IKaratzas and Shreve , 



1987)): 



V® = Ztt a (V t f + ! -^—(Ie(u) + /.(«)) du 

V Jtn &U-tQ 



(i) 

where Z\ is the process: 



Z[ l) = exp ( - (A< + t + gi o-iW t (i) ) . 



This models makes it natural to consider that the spike reception modifies the conductance of the network. 
This synapse model is considered first, and instantaneous current synapses next. 

A.2.1 Conductance-based synapses 

When neuron j receives a spike from one of its neighbors i, a current is generated, which has the value 
Wijgiyv) — V r ev) ■ Note that we artificially introduced V rev in the leak term, which amounts to formally 

(i) 

changing the current Ie , in order to integrate the equation in a simpler way. When neuron j receives a 
spike at time t* from neuron i, the conductance is hence increased by a coefficient Wij g. The solution of 
the membrane potential equation after time t* reads: 

V<*\t + t*) = VfrZ® + f 7«( S + OZE ds 
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where Z\ f = exp{— (Xj + -4- — Wijg)(t — t*) + CjW^ }. The reception of the spike therefore modifies 
the value of the membrane potential, which would have been equal to 

V^(t + t*) = VfcZte*"* 9 * + f li j) (s + t*)Z t - s e w » 9 ^- s) ds 

Jo 

if neuron j did not receive any spike. At time X*\ the value of (xi^ +t*) is equal to 9{X^ + t*), 
but after inhibition at time t* , the actual value of the membrane potential reads: 

V ij) {x¥ ] +f) =6e w v aX * + / I^{s + t*)Z t -s{e Wiias -l)ds 

Jo 

This value therefore depends on the history of the Brownian motion from time t* on. We are interested in 
finding the first spike time if nothing occurs in the network meanwhile. This is described by the random 
variable: 

THj = inf{< > t\ V t U) = 9{t)\V^(X { t P + t*) = 6(X ( t P + t*)} 

= X^i ] + mi{t > t*,V t U) = 0(t)\V u \X { t P +t*) = 0(X { t P + t*)} 



This case is therefore substantially more complex than the previous ones. Indeed, conditioning on the 
event that (t* + x[i') is the first crossing time of the threshold of conditions the trajectories of 
the Brownian motion in the interval [t* , t* + x[V], and we observe that the value of the updated process 
depends on the whole history of the Brownian motion in this interval (except in the very particular 
case where I e = 0, in which case the problem can be treated as previously). The problem of finding the 
law of the interaction random variable is therefore very tricky because of this complicated conditioning. 

A.2.2 Noisy Conductances and current-based synaptic interactions 

We now consider the case of instantaneous interactions. In that case, if neuron j receives a spike from 
neuron i at time t*, its membrane potential is instantaneously increased by a quantity equal to Wij. 
Similarly to the computations done in the previous case, the updated membrane potential now reads: 

v® (t + t*) = (V& + w tJ )z^ + f ip) (s + r)zE ds 

Jo 

and therefore at time X t . , this value reads: 

V«\xW+t*) = (6 + Wij )Z% 

(i) 

We therefore need to characterize the probability distribution of Z u) conditioned on the fact that the 

(i) 

neuron elicited a spike at time Xf, if it did not receive a spike at time t from neuron i. Here again, 
this problem depends in an intricate fashion on the distribution of the values of the underlying Brownian 
motion W^> at this time conditioned on the fact that the first crossing time of the process with the 
threshold was X y t t'. To the best of our knowledge, no solution is available so far. 
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A.3 Conclusion 

The current-based synapses are a necessary hypothesis in order to be able to express in closed-form the 
random variables involved in the transitions of the countdown process. The same problem arises in the 
case of nonlinear models: it remains very difficult to express the interaction random variable conditioned 
on the value of the countdown process at the time of the spike reception. 

It is important to note here that this limitation is also present in the deterministic case. Exten- 
sions of the framework of simple linear integrate-and-fire models are scarce and apply to very partic- 
ular models. The only successful extension to our knowledge was done by Tonnelier and collaborators 



dTonnelier et al 



20071) in the case of the quadratic integrate-and-fire neuron with constant input and 
constant threshold. It is based on the fact that the authors provide a closed-form expression for the 
membrane-potential. This fact is no more possible in the stochastic case, because even if we were to 
find a closed form expression for the membrane potential, the integration random variable will be very 
complex to express because of the conditioning on the value of the countdown variable as discussed 
previously. 

B Including Synaptic Delays and the Refractory Period 

In this appendix we include two biologically plausible phenomena in the description of the network 
activity, and see how this affects the Markovian framework. 

B.l Refractory period 

The refractory period is a transient phase just after the firing during which it is either impossible or 
difficult to communicate with the cell. This phenomenon is linked with the dynamics of ion channels 
and the hyperpolarization phase of the spike emission, lasts a few milliseconds, and prevents the neuron 
from firing spikes at an arbitrary high firing rate. It can be decomposed into two phases: the absolute 
refractory period, which is a constant period of time corresponding loosely to the hyperpolarization of 
the neuron during which is it impos sible to excite or inh ibit the cell no matter how great the stimulating 



current applied is (see for instance (|Ka ndel et al 



the phenomenon and (iGerstner and Kistlei 



2002a: 



2000, ch apter 9) for a further biological discussion of 



Arbib 



19981) for a discussion on the modeling this 



refractory period). Immediately after this phase begins the relative refractory period during which the 
initiation of a second action potential is inhibited but still possible. Modeling this relative refractory 
period amounts to considering that the synaptic inputs received at the level of the cell are weighted by a 
function depending on the time elapsed since the spike emission. 

For technical reasons, the relative refractory period is taken into account only for the spike inte- 
gration, and not for the noise integration. This assumption does not affect significantly the dynamics 
of the network, since the probability that the noisy current integrated during a time period as short as 
1 or 2ms be substantial is very small. For the spike integration this remark is no more valid, since a 
single spike induces large changes in the membrane voltage. During the relative refractory period, we 
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Figure 3: The refractory period at a spike emission, and the related k function weighting the synaptic 
inputs 
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Figure 4: There is a finite number of incoming spikes emitted before the postsynaptic cell fires. 



consider that synaptic efficiencies are weighted by a function depending on the time elapsed since the 
last spike was fired. We de note this function n{i), following the notation of Gerstner and Kistler in 
(IGerstner and Kistleru2002al) . In our case this function is unspecified, it is zero at t = and increases to 
1 with a characteristic time of around 2ms. It can be of bounded support or defined on M (see figure [3). 
We will see how taking in account this effect modifies the above instantaneous analysis, after introducing 
another essential effect occurring at the same time scale as the refractory period: synaptic delays. 



B.2 Synaptic delays 

Delay s are known t o be v ery important, for instance in shaping spatio-temporal dynamics of neuronal ac 



tivity dRoxin et al 



2005) or for generating global oscillations dBrunel and Hakim , 



1999). These synap- 



tic delays affect the interaction variable in a quite intricate fashion by adding a non-trivial memory-like 
phenomenon in the network. As a consequence of this we are not able to update instantaneously the 
countdown variable at each spike time: the transmission delays imply that one needs to keep the memory 
of a certain number of spikes times. Fortunately, because of the absolute refractory period, we only have 
to take into account a finite number of spikes that can possibly affect the postsynaptic potential after it 
elicits a spike (see figurelU). The maximal number of spikes concerned is given by M = maxjj Ll^J 
where Rj is the length of the absolute refractory period and [x\ is the floor function, i.e. the largest 
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integer smaller than or equal to x. 

Instead of considering the last firing times variables which contained the last firing time for each 
neuron of the network, we now consider the last M firing times variables. To this end, we define the 
matrix Hm G R WxM whose row i contains the M last firing times of neuron i. At the initial time, the 
M components of this row are set to the value miii,,{ //, — Ay}. If the neuron spikes at time £j, then 
each component of the row are modified: for all k e {2, . . . , N}, Hi^-i = and iJ^M = tj, and 
between two consecutive spikes of neuron i, the elements of the row i remain constant. 

Provided that the spike times are known, this matrix is a Markov chain. 

Indeed, let us denote by X n the countdown process of the network. The updates of the Markov chain 
(X, Hm) occur at spike emission times and spike arrivals at postsynaptic cells. The next spike if no 
delayed interaction occurs will be fired after a time given by r = min^ X™, and the first arrival of a 
possible spike at a cell is given by 

v = min {x = Hi k + An — t; . x > 0} 

<-.;• ! : v; ' 

fee{i, ...,M} 

• If this set is empty, the min is set to +00. If r < v, a spike is fired by the neuron i having the lowest 
countdown value. At this instant, the related countdown variable of neuron i is instantaneously 
reset by drawing in the law of the reset variable as described in section |2~T1 and row i of the matrix 
H is updated as indicated above. Other variables, such as the absolute time, are also to be updated 
at this point (e.g. the absolute time is updated to t m -\ + Xi, . ..). 

• If v < r, assume that the minimum is achieved for the value Hi,k+Aij for some i, j, k. This means 
that the k th latest spike of neuron i reaches the cell j and affects the cell in the same fashion as if 
there were no delay and a spike was emitted at this time. Therefore, the related interaction variable 
rjij of this connection will be added, and the countdown value of neuron j is updated together with 
the other Markov chain variables. The time is advanced to t + v. Note also that the min can be 
achieved at many different 3-tuples (i, j, k) at the same time. Moreover, it is also possible that 
an excitatory interaction makes a postsynaptic neuron fire instantaneously at the reception of the 
spike. All these cases can be treated sequentially, by iterating the mechanism we just described. 
Nevertheless, we are ensured that no avalanche can occur, because of the absolute refractory period 
and of the delays. 

B.3 The reset random variable 

Let us consider the effect of these features from the countdown process viewpoint. The reset variable is 
only affected by the absolute refractory period, and in a very simple way. Indeed, we formally consider 
that the neuron i is stuck at its reset value Vr during a fixed period of time Ri after having fired. After 
this period, the neuron membrane potential follows the evolution that depends on the particular model 
chosen. Therefore, the time of the next spike starting from time t + R4 has the same law as the reset 
variable in the case where we do not take into account the refractory period and the synaptic delay, i.e. 
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it has the law of the first hitting time, noted n, of the membrane potential process to the spike threshold 
with the time-shifted input I e (t + t* + Ri). The new reset variable of the related countdown process has 
simply the law of Yj = tj + Ri- 

The initialization random variable is clearly unchanged by taking into account synaptic delays, re- 
fractory periods and excitatory interactions. 

However, the case of the interaction variable is a little bit more intricate, as shown below. 

B.4 Interaction random variable 

Let us consider the effect of a spike emitted by neuron i and reaching cell j at time t in an inhibitory 
network. This spike will affect the membrane potential depending on the connectivity model chosen but 
corresponding to a synaptic weight Kj (t — tj)Wij where tj is the time of the last spike emitted by the cell 
J- 

Since we assumed the connection inhibitory, then the interaction random variable is readily deduced 
from the analysis of section|2]and of appendix[A] by changing the synaptic weight value ti>y to iuy Kj (ti — 
tj). Therefore in addition to the complementary variables necessary to define the countdown value, we 
need to keep in memory the last spike time for each neuron to define the synaptic weights and therefore 
the interaction random variable. Note that this random variable is almost surely equal to zero, whatever 
the neuron model considered, if the spike arrives during the absolute refractory period (i.e. t € [tj,tj + 
Rj]), since in that case n(t — tj) = 0. 

The spike transmission from the presynaptic cell to the postsynaptic one depends on the distance 
between the two cells, the speed of transmission of the signal along the axon and the transmission time at 
the synapse, and has a typical duration of a few milliseconds. To model the synaptic delay we consider 
that spikes emitted by a neuron affect the postsynaptic neurons after some delay Ay which can depend 
on both the presynaptic and the postsynaptic neurons. 

If neuron i fires a spike at time ti, its effect on the postsynaptic neuron j depends on the synaptic 
delay Ay, the countdown value of neuron j at this time X^\t{), and the time of the last spike emitted 
by j: 

(i) If Ay < X^'(t), then the reception of a spike at time U acts on the post-synaptic neuron at time 
t + Ay in the same way as discussed for the different models considered in sectionQ] but the effect 
can now be either excitatory or inhibitory, with a synaptic efficiency u>y Kj (ti + Ay —tj). 

(ii) If Ay > XV' (ti), the postsynaptic neuron will fire before receiving the spike from the presynaptic 
cell i, and it will act on the postsynaptic cell membrane with an efficiency WijKj(ti + Ay — Xj) 
at time ti + Ay . 
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